ICCB Environmental data download

Abstract

In this script we are downloading current and future environmental data (as rasters) to use in species distribution models (SDMs) for koalas (Phascolarctos cinereus) in the South-East Queensland (SEQ) region.

Import packages

library(terra)
library(dplyr)
library(sf)
library(ggplot2)

Load South East Queensland (SEQ) boundary

We start by defining our study area, which is the South East Queensland (SEQ) region. We will use the Local Government Areas (LGA) shapefile to define the extent of SEQ.

https://qldspatial.information.qld.gov.au/catalogue/custom/detail.page?fid={3F3DBD69-647B-4833-B0A5-CC43D5E70699}

# Load the study area shapefile
LGA <- st_read("Data/Environmental_variables/Local_Government_Areas.shp")

Reading layer `Local_Government_Areas' from data source 
  `/Users/scottforrest/Library/CloudStorage/OneDrive-QueenslandUniversityofTechnology/PhD - Scott Forrest/GIT/ICCB_geospatial_tools_conservation/Session 3/Data/Environmental_variables/Local_Government_Areas.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 78 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 137.9946 ymin: -29.17927 xmax: 153.5519 ymax: -9.087991
Geodetic CRS:  GDA94

# Check the coordinate reference system (CRS)
st_crs(LGA)

Coordinate Reference System:
  User input: GDA94 
  wkt:
GEOGCRS["GDA94",
    DATUM["Geocentric Datum of Australia 1994",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore."],
        BBOX[-60.55,93.41,-8.47,173.34]],
    ID["EPSG",4283]]

# Convert to WGS84
LGA <- LGA %>% st_transform(7856)

# Select local govt. areas for South East Queensland
LGA_SEQ <- LGA %>% 
  filter(lga %in% c("Brisbane City", 
                    "Moreton Bay City", 
                    "Logan City", 
                    "Ipswich City", 
                    "Redland City", 
                    "Scenic Rim Regional", 
                    "Somerset Regional", 
                    "Lockyer Valley Regional", 
                    "Gold Coast City", 
                    "Sunshine Coast Regional", 
                    "Toowoomba Regional", 
                    "Noosa Shire"))

ggplot() +
  geom_sf(data = LGA, color = "black") +
  geom_sf(data = LGA_SEQ, fill = "purple3", alpha = 0.5, color = "black", size = 0.2) +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(title = "Local Government Areas Queensland (SEQ in purple)")

ggplot() +
  geom_sf(data = LGA_SEQ, fill = "purple3", alpha = 0.5, color = "black", size = 0.2) +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(title = "Local Government Areas South East Queensland (SEQ)")

Merge into a single polygon

# Merge the SEQ LGAs into one polygon
SEQ_extent <- st_union(LGA_SEQ)

ggplot() +
  geom_sf(data = SEQ_extent, fill = "purple3", alpha = 0.5, color = "black", size = 0.2) +
  theme_minimal() +
  theme(legend.position = "none") +
  ggtitle("South-East Queensland Spatial Extent") + 
  theme_bw() 

Save SEQ extent for other scripts

# Convert our SEQ extent to a SpatExtent object by converting to a SpatVector
SEQ_extent.vect <- terra::vect(SEQ_extent)

# Write the SEQ extent to a shapefile
writeVector(SEQ_extent.vect, "Data/Environmental_variables/SEQ_extent.shp", overwrite = T)

Load current environmental data

Layers were made available to us by the EcoCommons team and were created by Toombs and Ma (2025):

Toombs, N., and Ma S., 2025, A High-Resolution Dataset of 19 Bioclimatic Indices over Australia, Climate Projections and Services – Queensland Treasury, Brisbane, Queensland. [https://longpaddock.qld.gov.au/qld-future-climate/data-info/tern/]

files <- list.files("Data/Environmental_variables/Current_climate_QLD", 
             pattern = ".tif$", 
             full.names = TRUE)

# Load all bioclim rasters
current_bioclim <- lapply(files, terra::rast) 

# Make into one raster stack
current_bioclim <- rast(current_bioclim)

# Plot the current bioclimatic variables
plot(current_bioclim)

# Examine the resolution
current_bioclim

class       : SpatRaster 
dimensions  : 681, 841, 19  (nrow, ncol, nlyr)
resolution  : 0.05, 0.05  (x, y)
extent      : 111.975, 154.025, -44.025, -9.975  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : bioclim_01.tif  
              bioclim_02.tif  
              bioclim_03.tif  
              ... and 16 more sources
names       : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ... 

# Check the CRS
crs(current_bioclim)

[1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        MEMBER[\"World Geodetic System 1984 (G2296)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

# Update CRS to EPSG:7856 (GDA2020 / MGA zone 56)
current_bioclim <- terra::project(current_bioclim, "EPSG:7856")

# Our resolution is now ~5km by 5km
current_bioclim

class       : SpatRaster 
dimensions  : 834, 898, 19  (nrow, ncol, nlyr)
resolution  : 5590.925, 5590.925  (x, y)
extent      : -4407995, 612655.7, 4234346, 8897178  (xmin, xmax, ymin, ymax)
coord. ref. : GDA2020 / MGA zone 56 (EPSG:7856) 
source(s)   : memory
names       : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ... 
min values  :   5.182303,   5.520304,   35.25507,   105.7126,   16.82971,  -4.954343, ... 
max values  :  29.453394,  16.870607,   61.84691,   680.0889,   42.45033,  22.998774, ... 

Mask and crop to SEQ extent

# Mask the current bioclimatic variables to the SEQ extent
current_bioclim <- terra::mask(current_bioclim, SEQ_extent.vect)

# Crop to the SEQ region
current_bioclim <- terra::crop(current_bioclim, SEQ_extent.vect)

# Plot all current bioclimatic variables
plot(current_bioclim)

# Save the current environmental covariates
writeRaster(current_bioclim,
            filename = "Data/Environmental_variables/SEQ_current_bioclim.tif",
            overwrite = T)

Load future environmental data

Here we load outputs from a moderate-high emissions shared socio-economic path scenario (SSP 3.70) for the year 2090 (2080 - 2099).

files <- list.files("Data/Environmental_variables/Future_climate_SSP370_2090", 
             pattern = ".tif$", 
             full.names = TRUE)

# Load all bioclim rasters
future_bioclim <- lapply(files, terra::rast) 

# Make into one raster stack
future_bioclim <- rast(future_bioclim)

# Plot the future bioclimatic variables
plot(future_bioclim)

# Examine the resolution
future_bioclim

class       : SpatRaster 
dimensions  : 681, 841, 19  (nrow, ncol, nlyr)
resolution  : 0.05, 0.05  (x, y)
extent      : 111.975, 154.025, -44.025, -9.975  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : bioclim_01.tif  
              bioclim_02.tif  
              bioclim_03.tif  
              ... and 16 more sources
names       : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ... 

# Check the CRS
crs(future_bioclim)

[1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        MEMBER[\"World Geodetic System 1984 (G2296)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

# Update CRS 
future_bioclim <- terra::project(future_bioclim, "EPSG:7856")

# Our resolution is now ~5km by 5km
future_bioclim

class       : SpatRaster 
dimensions  : 834, 898, 19  (nrow, ncol, nlyr)
resolution  : 5590.925, 5590.925  (x, y)
extent      : -4407995, 612655.7, 4234346, 8897178  (xmin, xmax, ymin, ymax)
coord. ref. : GDA2020 / MGA zone 56 (EPSG:7856) 
source(s)   : memory
names       : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ... 
min values  :    0.00000,    0.00000,     0.0000,     0.0000,    0.00000,  -2.243844, ... 
max values  :   33.02249,   16.09955,    66.1461,   706.1523,   46.24717,  25.827158, ... 

Mask to SEQ extent

We can see that for these layers the water surrounding Australia is not NA, but are values of 0. In some cells this introduces an artifact on the coastline where the values in some cells are an average between realistic values and 0, which results in cells with unusual and unrealistic values. We will do our best to mask these out.

# First mask and crop by the SEQ extent
future_bioclim <- terra::mask(future_bioclim, SEQ_extent.vect) 
future_bioclim <- terra::crop(future_bioclim, SEQ_extent.vect)

# Plot one of the variables - max temp of the warmest month
plot(future_bioclim[[5]], main = "Future BIO5")

We can clearly see the artifacts in BIO5, which are not present in the current layers.

plot(current_bioclim[[5]], main = "Current BIO5")

If we plot a histogram of the values we can see where we can set a threshold to exclude those lower values on the coastline. We tested values to ensure that we were removing as many of the dodgy values at the coast whilst retaining all real values, and the best balance was 28 degrees.

Unfortunately now we have slightly fewer cells on the coastline, but as these were artifacts we don’t have the true data for those cells anyway.

# plot the distribution of values in BIO5
hist(future_bioclim[[5]], breaks = 100)

# create a mask for the artifacts, which we will apply across all of the layers
artifact_mask <- future_bioclim[[5]] > 28
names(artifact_mask) <- "artifact_mask"

# plot the artifact mask
plot(artifact_mask)

# plot the updated future BIO5 layer
plot(future_bioclim[[5]], main = "Future BIO5 - unmasked")

# set all artifact values to NA (across all layers)
future_bioclim <- terra::mask(future_bioclim, artifact_mask, maskvalues = 0) # Set all values of 0 to NA

# plot the updated future BIO5 layer
plot(future_bioclim[[5]], main = "Future BIO5 - masked")

plot(current_bioclim[[5]], main = "Current BIO5")

# this removes those lower values
hist(future_bioclim[[5]], breaks = 100)

We can see there are a couple raster cells that have artifacts, but we can filter by the difference between the current and future layers to remove these. Unfortunately, temperature is only going to increase, so we can filter by values where the difference is negative for the BIOCLIM 5 layer, which is the max temperature of the warmest month.

# have a look at the difference between current and future
plot(future_bioclim[[5]] - current_bioclim[[5]])

min(values(future_bioclim[[5]] - current_bioclim[[5]]), na.rm = TRUE)

[1] -1.306969

# save that and create a mask for the artifact cells
temp_diff <- future_bioclim[[5]] - current_bioclim[[5]]
temp_diff_mask <- temp_diff > 1 # there is one remaining artifact cell that's about 0.7, so we'll remove that

# have a look at the mask
plot(temp_diff_mask)

# these values look much better without the artifacts
plot(mask(future_bioclim[[5]], temp_diff_mask, maskvalues = 0), main = "Future BIO5 - masked by difference")

# we'll apply that mask to all layers
future_bioclim <- mask(future_bioclim, temp_diff_mask, maskvalues = 0)

Now that we have applied the artifact mask across all layers, they look better.

plot(future_bioclim)

Now we can save the future environmental covariates for SEQ.

# Save the cropped and masked future environmental covariates
writeRaster(future_bioclim,
            filename = "Data/Environmental_variables/SEQ_future_bioclim.2090.SSP370.tif",
            overwrite = T)

Load future environmental data 2

Here we load outputs from a low emissions shared socio-economic path scenario (SSP 1.26) for the year 2090 (2080 - 2099).

files <- list.files("Data/Environmental_variables/Future_climate_SSP126_2090", 
             pattern = ".tif$", 
             full.names = TRUE)

# Load all bioclim rasters
future_bioclim <- lapply(files, terra::rast) 

# Make into one raster stack
future_bioclim <- rast(future_bioclim)

plot(future_bioclim)

# Examine the resolution
future_bioclim

class       : SpatRaster 
dimensions  : 681, 841, 19  (nrow, ncol, nlyr)
resolution  : 0.05, 0.05  (x, y)
extent      : 111.975, 154.025, -44.025, -9.975  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : bioclim_01.tif  
              bioclim_02.tif  
              bioclim_03.tif  
              ... and 16 more sources
names       : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ... 
min values  :    0.00000,    0.00000,     0.0000,         ? ,         ? ,         ? , ... 
max values  :   30.47033,   16.75114,    63.2751,         ? ,         ? ,         ? , ... 

# Check the CRS
crs(future_bioclim)

[1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        MEMBER[\"World Geodetic System 1984 (G2296)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

# Update CRS 
future_bioclim <- terra::project(future_bioclim, "EPSG:7856")

# Our resolution is now ~5km by 5km
future_bioclim

class       : SpatRaster 
dimensions  : 834, 898, 19  (nrow, ncol, nlyr)
resolution  : 5590.925, 5590.925  (x, y)
extent      : -4407995, 612655.7, 4234346, 8897178  (xmin, xmax, ymin, ymax)
coord. ref. : GDA2020 / MGA zone 56 (EPSG:7856) 
source(s)   : memory
names       : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ... 
min values  :    0.00000,    0.00000,    0.00000,     0.0000,    0.00000,  -4.081378, ... 
max values  :   30.78085,   16.89395,   63.34254,   687.5474,   43.83962,  23.968334, ... 

Mask to SEQ extent

We will crop and mask the future bioclimatic variables to the SEQ extent, and apply the artifact mask as we did for the SSP 3.70 layers.

# First mask and crop by the SEQ extent
future_bioclim <- terra::mask(future_bioclim, SEQ_extent.vect) 
future_bioclim <- terra::crop(future_bioclim, SEQ_extent.vect)

# set all artifact values to NA
future_bioclim <- terra::mask(future_bioclim, artifact_mask, maskvalues = 0) # Set all values of 0 to NA
future_bioclim <- terra::mask(future_bioclim, temp_diff_mask, maskvalues = 0) # also use the temp difference mask
plot(future_bioclim[[5]])

plot(future_bioclim)

# Save the future environmental covariates
writeRaster(future_bioclim,
            filename = "Data/Environmental_variables/SEQ_future_bioclim.2090.SSP126.tif",
            overwrite = T)